Dbike    = (CBD_Data$MainModeENG == "Bicycle")
Dwalk    = (CBD_Data$MainModeENG == "Walk")
DPT      = (CBD_Data$MainModeENG == "Rail" | CBD_Data$MainModeENG == "Bus")
DNM      = (CBD_Data$MainModeENG == "Bicycle" | CBD_Data$MainModeENG == "Walk")
DIntra   = ((CBD_Data$DistrictCodeO == 103&CBD_Data$DistrictCodeD == 103)|(CBD_Data$DistrictCodeO == 104&CBD_Data$DistrictCodeD == 104))
DInter   = !DIntra
## Calculate the Log-likelihood function
LL = colSums(Dtrain*log(Ptrain1) + Dbus*log(Pbus1) + Dcar*log(Pcar1) + Dbike*log(Pbike1) + Dwalk*log(Pwalk1) +
Dtrain*log(Ptrain2) + Dbus*log(Pbus2) + Dcar*log(Pcar2) + Dbike*log(Pbike2) + Dwalk*log(Pwalk2)
# DPT*log(P.PT1) + DNM*log(P.NM1) +
# DPT*log(P.PT2) + DNM*log(P.NM2) +
# DIntra*log(P.Intra) + DInter*log(P.Inter)
)
}
##### Maximize the Log-likelihood function#####
## Setting up initValaue
NL_initValue = list(Ctrain = 0,
Cbus   = 0,
Ccar   = 0,
Cbike  = 0,
TTime  = 0,
TCost  = 0,
WDayTrain = 0,
WDayBus = 0,
WDayCar = 0,
WDayBike = 0,
WDayWalk = 0,
WPurTrain = 0,
WPurBus = 0,
WPurCar = 0,
WPurBike = 0,
WPurWalk = 0,
LogsumNM = 1,
LogsumPT = 1,
LogsumInter = 1,
LogsumIntra = 1
# # Logsum1 = 1
# Logsum2 = 1
)
# Clear the Console in R
cat("\014")
## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue, method = "BFGS", control = list(trace = 1000000, maxit = 100000000, fnscale = -1))
NL_LL0 = N*log(1/5)     ## Initial log-likelihood
NL_LL1 = logLik(NL)             ## Converged log-likelihood
## McFadden's Rho square
NL_Rho     = (1 - (NL_LL1/NL_LL0))
NL_Rho.adj = (1 - ((NL_LL1 - length(coef(NL)))/NL_LL0))
NL_GoF = cbind(NL_LL0,NL_LL1, NL_Rho, NL_Rho.adj)  ## Goodness of fit
## Preview NL model result in console
print(summary(NL))
print(NL_GoF)
### Multinomial Logit model estimation
rm(list = ls())
source("BIN2016-TeamC-MNL.R")
##### Define the log-likelihood function of the logit model#####
NL_LL <- function(Ctrain,
Cbus,
Ccar,
Cbike,
TTime,
TCost,
WDayTrain,
WDayBus,
WDayCar,
WDayBike,
WDayWalk,
WPurTrain,
WPurBus,
WPurCar,
WPurBike,
WPurWalk,
LogsumPT,
LogsumNM,
LogsumIntra,
LogsumInter
# # Logsum1,
# Logsum2
)
{
# LogsumPT = Logsum1
# LogsumNM = Logsum1
# LogsumIntra = Logsum2
# LogsumInter = Logsum2
## Declare the log-likelihood variable, set value to 0
NL_LL = 0
## Utility functions
Vtrain  = Ctrain*matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeTrain  + TCost *CBD_Data$FareTrain                        + WDayTrain *CBD_Data$DWeekday + WPurTrain *CBD_Data$DWPur
Vbus    = Cbus  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeBus    + TCost *CBD_Data$FareBus                          + WDayBus *CBD_Data$DWeekday   + WPurBus *CBD_Data$DWPur
Vcar    = Ccar  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeCar         + TCost *(CBD_Data$ParkingCost + CBD_Data$FareCar) + WDayCar *CBD_Data$DWeekday   + WPurCar *CBD_Data$DWPur
Vbike   = Cbike *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeBike                                                           + WDayBike *CBD_Data$DWeekday  + WPurBike *CBD_Data$DWPur
Vwalk   =        matrix(0,nrow = N, ncol = 1) + TTime *CBD_Data$TimeWalk                                                           + WDayWalk *CBD_Data$DWeekday  + WPurWalk *CBD_Data$DWPur
V.PT1   = 2
V.NM1   = 2
V.PT2   = 2
V.NM2   = 2
V.Intra = 2
V.Inter = 2
Vtrain = CBD_Data$ModeAvailableTrain * Vtrain
Vbus   = CBD_Data$ModeAvailableBus * Vbus
Vcar   = CBD_Data$ModeAvailableCar * Vcar
Vbike  = CBD_Data$ModeAvailableBike * Vbike
Vwalk  = CBD_Data$ModeAvailableWalk * Vwalk
## Log-sum variables
LS.PT1 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
LS.NM1 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
LS.PT2 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
LS.NM2 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
LS.Intra = log(exp(LS.PT1/LogsumIntra) + exp(LS.NM1/LogsumIntra))
LS.Inter = log(exp(LS.PT2/LogsumInter) + exp(LS.NM2/LogsumInter))
## Denominators
DenoPT1  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
DenoNM1  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)
DenoPT2  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
DenoNM2  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)
DenoIntra = exp(Vcar) + exp(V.PT1 + LogsumPT*LS.PT1) + exp(V.NM1 + LogsumNM*LS.NM1)
DenoInter = exp(Vcar) + exp(V.PT2 + LogsumPT*LS.PT2) + exp(V.NM2 + LogsumNM*LS.NM2)
Deno = exp(V.Intra + LogsumIntra*LS.Intra) + exp(V.Inter + LogsumInter*LS.Inter)
## Probabilites
# Conditional probabilits
Ptrain.PT1  = exp(Vtrain/LogsumPT)/DenoPT1
Pbus.PT1    = exp(Vbus/LogsumPT)/DenoPT1
Pbike.NM1   = exp(Vbike/LogsumNM)/DenoNM1
Pwalk.NM1   = exp(Vwalk/LogsumNM)/DenoNM1
Ptrain.PT2  = exp(Vtrain/LogsumPT)/DenoPT2
Pbus.PT2    = exp(Vbus/LogsumPT)/DenoPT2
Pbike.NM2   = exp(Vbike/LogsumNM)/DenoNM2
Pwalk.NM2   = exp(Vwalk/LogsumNM)/DenoNM2
Pcar.Intra  = exp(Vcar)/DenoIntra
P.PT1       = exp(V.PT1 + LogsumPT*LS.PT1)/DenoIntra
P.NM1       = exp(V.NM1 + LogsumNM*LS.NM1)/DenoIntra
Pcar.Inter  = exp(Vcar)/DenoInter
P.PT2       = exp(V.PT2 + LogsumPT*LS.PT2)/DenoInter
P.NM2       = exp(V.NM2 + LogsumNM*LS.NM2)/DenoInter
P.Intra     = exp(V.Intra + LogsumIntra*LS.Intra)/Deno
P.Inter     = exp(V.Inter + LogsumInter*LS.Inter)/Deno
##Joint probabilities
Ptrain1     = Ptrain.PT1*P.PT1*P.Intra
Pbus1       = Pbus.PT1*P.PT1*P.Intra
Pbike1      = Pbike.NM1*P.NM1*P.Intra
Pwalk1      = Pwalk.NM1*P.NM1*P.Intra
Pcar1       = Pcar.Intra*P.Intra
Ptrain2     = Ptrain.PT2*P.PT2*P.Inter
Pbus2       = Pbus.PT2*P.PT2*P.Inter
Pbike2      = Pbike.NM2*P.NM2*P.Inter
Pwalk2      = Pwalk.NM2*P.NM2*P.Inter
Pcar2       = Pcar.Inter*P.Inter
Ptrain1 = (Ptrain1 != 0)*Ptrain1 + (Ptrain1 == 0)
Pbus1   = (Pbus1   != 0)*Pbus1   + (Pbus1   == 0)
Pcar1   = (Pcar1   != 0)*Pcar1   + (Pcar1   == 0)
Pbike1  = (Pbike1  != 0)*Pbike1  + (Pbike1  == 0)
Pwalk1  = (Pwalk1  != 0)*Pwalk1  + (Pwalk1  == 0)
Ptrain2 = (Ptrain2 != 0)*Ptrain2 + (Ptrain2 == 0)
Pbus2   = (Pbus2   != 0)*Pbus2   + (Pbus2   == 0)
Pcar2   = (Pcar2   != 0)*Pcar2   + (Pcar2   == 0)
Pbike2  = (Pbike2  != 0)*Pbike2  + (Pbike2  == 0)
Pwalk2  = (Pwalk2  != 0)*Pwalk2  + (Pwalk2  == 0)
## Dummies for choice
Dtrain   = (CBD_Data$MainModeENG == "Rail")
Dbus     = (CBD_Data$MainModeENG == "Bus")
Dcar     = (CBD_Data$MainModeENG == "Car")
Dbike    = (CBD_Data$MainModeENG == "Bicycle")
Dwalk    = (CBD_Data$MainModeENG == "Walk")
DPT      = (CBD_Data$MainModeENG == "Rail" | CBD_Data$MainModeENG == "Bus")
DNM      = (CBD_Data$MainModeENG == "Bicycle" | CBD_Data$MainModeENG == "Walk")
DIntra   = ((CBD_Data$DistrictCodeO == 103&CBD_Data$DistrictCodeD == 103)|(CBD_Data$DistrictCodeO == 104&CBD_Data$DistrictCodeD == 104))
DInter   = !DIntra
## Calculate the Log-likelihood function
LL = -colSums(Dtrain*log(Ptrain1) + Dbus*log(Pbus1) + Dcar*log(Pcar1) + Dbike*log(Pbike1) + Dwalk*log(Pwalk1) +
Dtrain*log(Ptrain2) + Dbus*log(Pbus2) + Dcar*log(Pcar2) + Dbike*log(Pbike2) + Dwalk*log(Pwalk2)
# DPT*log(P.PT1) + DNM*log(P.NM1) +
# DPT*log(P.PT2) + DNM*log(P.NM2) +
# DIntra*log(P.Intra) + DInter*log(P.Inter)
)
}
##### Maximize the Log-likelihood function#####
## Setting up initValaue
NL_initValue = list(Ctrain = 0,
Cbus   = 0,
Ccar   = 0,
Cbike  = 0,
TTime  = 0,
TCost  = 0,
WDayTrain = 0,
WDayBus = 0,
WDayCar = 0,
WDayBike = 0,
WDayWalk = 0,
WPurTrain = 0,
WPurBus = 0,
WPurCar = 0,
WPurBike = 0,
WPurWalk = 0,
LogsumNM = 1,
LogsumPT = 1,
LogsumInter = 1,
LogsumIntra = 1
# # Logsum1 = 1
# Logsum2 = 1
)
# Clear the Console in R
cat("\014")
## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue, method = "BFGS", control = list(trace = 1000000, maxit = 100000000))
NL_LL0 = N*log(1/5)     ## Initial log-likelihood
NL_LL1 = logLik(NL)             ## Converged log-likelihood
## McFadden's Rho square
NL_Rho     = (1 - (NL_LL1/NL_LL0))
NL_Rho.adj = (1 - ((NL_LL1 - length(coef(NL)))/NL_LL0))
NL_GoF = cbind(NL_LL0,NL_LL1, NL_Rho, NL_Rho.adj)  ## Goodness of fit
## Preview NL model result in console
print(summary(NL))
print(NL_GoF)
NL = mle2(NL_LL, start = NL_initValue,  control = list(trace = 1000000, maxit = 100000000))
NL = mle2(NL_LL, start = NL_initValue,  control = list(trace = F, maxit = 100000000))
NL = mle2(NL_LL, start = NL_initValue,  control = list(trace = T, maxit = 100000000))
source("BIN2016-TeamC-MNL.R")
source("BIN2016-TeamC-MNL.R")
######################################################################
#                                                                    #
# Name  :  15th Travel Behavior and Network Model Course             #
# Author:  Team C                                                    #
#          Canh Do, Hong Nguyen, Aktar Hossain, Perez Barbosa David  #
#          Transportation Study Group (TSG),Hiroshima University     #
# Time  :  23th-25th September, 2016                                 #
#                                                                    #
######################################################################
#Run MNL Code
rm(list = ls())
source("BIN2016-TeamC-MNL.R")
##### Define the log-likelihood function of the logit model#####
NL_LL <- function(Ctrain,
Cbus,
Ccar,
Cbike,
TTime,
TCost,
WDayTrain,
WDayBus,
WDayCar,
WDayBike,
WDayWalk,
WPurTrain,
WPurBus,
WPurCar,
WPurBike,
WPurWalk,
LogsumPT,
LogsumNM
)
{
## Declare the log-likelihood variable, set value to 0
NL_LL = 0
## Utility functions
Vtrain  = Ctrain*matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeTrain  + TCost *CBD_Data$FareTrain                        + WDayTrain *CBD_Data$DWeekday + WPurTrain *CBD_Data$DWPur
Vbus    = Cbus  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeBus    + TCost *CBD_Data$FareBus                          + WDayBus *CBD_Data$DWeekday   + WPurBus *CBD_Data$DWPur
Vcar    = Ccar  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeCar         + TCost *(CBD_Data$ParkingCost + CBD_Data$FareCar) + WDayCar *CBD_Data$DWeekday   + WPurCar *CBD_Data$DWPur
Vbike   = Cbike *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeBike                                                           + WDayBike *CBD_Data$DWeekday  + WPurBike *CBD_Data$DWPur
Vwalk   =        matrix(0,nrow = N, ncol = 1) + TTime *CBD_Data$TimeWalk                                                           + WDayWalk *CBD_Data$DWeekday  + WPurWalk *CBD_Data$DWPur
V.PT    = 2
V.NM    = 2
Vtrain = CBD_Data$ModeAvailableTrain * Vtrain
Vbus   = CBD_Data$ModeAvailableBus * Vbus
Vcar   = CBD_Data$ModeAvailableCar * Vcar
Vbike  = CBD_Data$ModeAvailableBike * Vbike
Vwalk  = CBD_Data$ModeAvailableWalk * Vwalk
## Log-sum variables
LS.PT = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
LS.NM = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
## Denominators
DenoPT  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
DenoNM  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)
Deno = exp(Vcar) + exp(V.PT + LogsumPT*LS.PT) + exp(V.NM + LogsumNM*LS.NM)
## Probabilites
# Conditional probabilits
Ptrain.PT  = exp(Vtrain/LogsumPT)/DenoPT
Pbus.PT    = exp(Vbus  /LogsumPT)/DenoPT
Pbike.NM   = exp(Vbike /LogsumNM)/DenoNM
Pwalk.NM   = exp(Vwalk /LogsumNM)/DenoNM
P.PT       = exp(V.PT + LogsumPT*LS.PT)/Deno
P.NM       = exp(V.NM + LogsumNM*LS.NM)/Deno
Pcar       = exp(Vcar)/Deno
##Joint probabilities
Ptrain     = Ptrain.PT*P.PT
Pbus       = Pbus.PT*P.PT
Pbike      = Pbike.NM*P.NM
Pwalk      = Pwalk.NM*P.NM
## Avoid problems stemming from choice probabilities becoming zero.
Ptrain = (Ptrain != 0)*Ptrain + (Ptrain == 0)
Pbus   = (Pbus   != 0)*Pbus   + (Pbus   == 0)
Pcar   = (Pcar   != 0)*Pcar   + (Pcar   == 0)
Pbike  = (Pbike  != 0)*Pbike  + (Pbike  == 0)
Pwalk  = (Pwalk  != 0)*Pwalk  + (Pwalk  == 0)
## Dummies for mode choice
Dtrain   = (CBD_Data$MainModeENG == "Rail")
Dbus     = (CBD_Data$MainModeENG == "Bus")
Dcar     = (CBD_Data$MainModeENG == "Car")
Dbike    = (CBD_Data$MainModeENG == "Bicycle")
Dwalk    = (CBD_Data$MainModeENG == "Walk")
DPT      = (CBD_Data$MainModeENG == "Rail" | CBD_Data$MainModeENG == "Bus")
DNM      = (CBD_Data$MainModeENG == "Bicycle" | CBD_Data$MainModeENG == "Walk")
## Calculate the Log-likelihood function
NL_LL = -colSums(Dtrain*log(Ptrain) + Dbus*log(Pbus) + Dcar*log(Pcar) + Dbike*log(Pbike) + Dwalk*log(Pwalk) #+
#DPT*log(P.PT) + DNM*log(P.NM)
)
return(NL_LL)
}
##### Maximize the Log-likelihood function#####
## Setting up initValaue
NL_initValue = list(Ctrain = 0,
Cbus   = 0,
Ccar   = 0,
Cbike  = 0,
TTime  = 0,
TCost  = 0,
WDayTrain = 0,
WDayBus = 0,
WDayCar = 0,
WDayBike = 0,
WDayWalk = 0,
WPurTrain = 0,
WPurBus = 0,
WPurCar = 0,
WPurBike = 0,
WPurWalk = 0,
LogsumPT = 1,
LogsumNM = 1)
## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue, method = "BFGS", control = list(trace = T, maxit = 100000000))
NL_LL0 = N*log(1/5)     ## Initial log-likelihood
NL_LL1 = logLik(NL)             ## Converged log-likelihood
## McFadden's Rho square
NL_Rho     = (1 - (NL_LL1/NL_LL0))
NL_Rho.adj = (1 - ((NL_LL1 - length(coef(NL)))/NL_LL0))
NL_GoF = cbind(NL_LL0,NL_LL1, NL_Rho, NL_Rho.adj)  ## Goodness of fit
## Preview NL model result in console
print(summary(NL))
print(NL_GoF)
######################################################################
#                                                                    #
# Name  :  15th Travel Behavior and Network Model Course             #
# Author:  Team C                                                    #
#          (Canh Do, Hong Nguyen, Aktar Hossain, Perez Barbosa David)#
#          Transportation Study Group (TSG),Hiroshima University     #
# Time  :  23th-25th September, 2016                                 #
#                                                                    #
######################################################################
######################################################################
#Run MNL Code
rm(list = ls())
source("BIN2016-TeamC-MNL.R")
##### Define the log-likelihood function of the logit model#####
NL_LL <- function(Ctrain,
Cbus,
Ccar,
Cbike,
TTime,
TCost,
WDayTrain,
WDayBus,
WDayCar,
WDayBike,
WDayWalk,
WPurTrain,
WPurBus,
WPurCar,
WPurBike,
WPurWalk,
LogsumPT,
LogsumNM,
LogsumIntra,
LogsumInter
# # Logsum1,
# Logsum2
)
{
# LogsumPT = Logsum1
# LogsumNM = Logsum1
# LogsumIntra = Logsum2
# LogsumInter = Logsum2
## Declare the log-likelihood variable, set value to 0
NL_LL = 0
## Utility functions
Vtrain  = Ctrain*matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeTrain  + TCost *CBD_Data$FareTrain                        + WDayTrain *CBD_Data$DWeekday + WPurTrain *CBD_Data$DWPur
Vbus    = Cbus  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TotalTimeBus    + TCost *CBD_Data$FareBus                          + WDayBus *CBD_Data$DWeekday   + WPurBus *CBD_Data$DWPur
Vcar    = Ccar  *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeCar         + TCost *(CBD_Data$ParkingCost + CBD_Data$FareCar) + WDayCar *CBD_Data$DWeekday   + WPurCar *CBD_Data$DWPur
Vbike   = Cbike *matrix(1,nrow = N, ncol = 1) + TTime *CBD_Data$TimeBike                                                           + WDayBike *CBD_Data$DWeekday  + WPurBike *CBD_Data$DWPur
Vwalk   =        matrix(0,nrow = N, ncol = 1) + TTime *CBD_Data$TimeWalk                                                           + WDayWalk *CBD_Data$DWeekday  + WPurWalk *CBD_Data$DWPur
V.PT1   = 2
V.NM1   = 2
V.PT2   = 2
V.NM2   = 2
V.Intra = 2
V.Inter = 2
Vtrain = CBD_Data$ModeAvailableTrain * Vtrain
Vbus   = CBD_Data$ModeAvailableBus * Vbus
Vcar   = CBD_Data$ModeAvailableCar * Vcar
Vbike  = CBD_Data$ModeAvailableBike * Vbike
Vwalk  = CBD_Data$ModeAvailableWalk * Vwalk
## Log-sum variables
LS.PT1 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
LS.NM1 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
LS.PT2 = log(exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT))
LS.NM2 = log(exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM))
LS.Intra = log(exp(LS.PT1/LogsumIntra) + exp(LS.NM1/LogsumIntra))
LS.Inter = log(exp(LS.PT2/LogsumInter) + exp(LS.NM2/LogsumInter))
## Denominators
DenoPT1  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
DenoNM1  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)
DenoPT2  = exp(Vtrain/LogsumPT) + exp(Vbus/LogsumPT)
DenoNM2  = exp(Vbike/LogsumNM)  + exp(Vwalk/LogsumNM)
DenoIntra = exp(Vcar) + exp(V.PT1 + LogsumPT*LS.PT1) + exp(V.NM1 + LogsumNM*LS.NM1)
DenoInter = exp(Vcar) + exp(V.PT2 + LogsumPT*LS.PT2) + exp(V.NM2 + LogsumNM*LS.NM2)
Deno = exp(V.Intra + LogsumIntra*LS.Intra) + exp(V.Inter + LogsumInter*LS.Inter)
## Probabilites
# Conditional probabilits
Ptrain.PT1  = exp(Vtrain/LogsumPT)/DenoPT1
Pbus.PT1    = exp(Vbus/LogsumPT)/DenoPT1
Pbike.NM1   = exp(Vbike/LogsumNM)/DenoNM1
Pwalk.NM1   = exp(Vwalk/LogsumNM)/DenoNM1
Ptrain.PT2  = exp(Vtrain/LogsumPT)/DenoPT2
Pbus.PT2    = exp(Vbus/LogsumPT)/DenoPT2
Pbike.NM2   = exp(Vbike/LogsumNM)/DenoNM2
Pwalk.NM2   = exp(Vwalk/LogsumNM)/DenoNM2
Pcar.Intra  = exp(Vcar)/DenoIntra
P.PT1       = exp(V.PT1 + LogsumPT*LS.PT1)/DenoIntra
P.NM1       = exp(V.NM1 + LogsumNM*LS.NM1)/DenoIntra
Pcar.Inter  = exp(Vcar)/DenoInter
P.PT2       = exp(V.PT2 + LogsumPT*LS.PT2)/DenoInter
P.NM2       = exp(V.NM2 + LogsumNM*LS.NM2)/DenoInter
P.Intra     = exp(V.Intra + LogsumIntra*LS.Intra)/Deno
P.Inter     = exp(V.Inter + LogsumInter*LS.Inter)/Deno
##Joint probabilities
Ptrain1     = Ptrain.PT1*P.PT1*P.Intra
Pbus1       = Pbus.PT1*P.PT1*P.Intra
Pbike1      = Pbike.NM1*P.NM1*P.Intra
Pwalk1      = Pwalk.NM1*P.NM1*P.Intra
Pcar1       = Pcar.Intra*P.Intra
Ptrain2     = Ptrain.PT2*P.PT2*P.Inter
Pbus2       = Pbus.PT2*P.PT2*P.Inter
Pbike2      = Pbike.NM2*P.NM2*P.Inter
Pwalk2      = Pwalk.NM2*P.NM2*P.Inter
Pcar2       = Pcar.Inter*P.Inter
Ptrain1 = (Ptrain1 != 0)*Ptrain1 + (Ptrain1 == 0)
Pbus1   = (Pbus1   != 0)*Pbus1   + (Pbus1   == 0)
Pcar1   = (Pcar1   != 0)*Pcar1   + (Pcar1   == 0)
Pbike1  = (Pbike1  != 0)*Pbike1  + (Pbike1  == 0)
Pwalk1  = (Pwalk1  != 0)*Pwalk1  + (Pwalk1  == 0)
Ptrain2 = (Ptrain2 != 0)*Ptrain2 + (Ptrain2 == 0)
Pbus2   = (Pbus2   != 0)*Pbus2   + (Pbus2   == 0)
Pcar2   = (Pcar2   != 0)*Pcar2   + (Pcar2   == 0)
Pbike2  = (Pbike2  != 0)*Pbike2  + (Pbike2  == 0)
Pwalk2  = (Pwalk2  != 0)*Pwalk2  + (Pwalk2  == 0)
## Dummies for choice
Dtrain   = (CBD_Data$MainModeENG == "Rail")
Dbus     = (CBD_Data$MainModeENG == "Bus")
Dcar     = (CBD_Data$MainModeENG == "Car")
Dbike    = (CBD_Data$MainModeENG == "Bicycle")
Dwalk    = (CBD_Data$MainModeENG == "Walk")
DPT      = (CBD_Data$MainModeENG == "Rail" | CBD_Data$MainModeENG == "Bus")
DNM      = (CBD_Data$MainModeENG == "Bicycle" | CBD_Data$MainModeENG == "Walk")
DIntra   = ((CBD_Data$DistrictCodeO == 103&CBD_Data$DistrictCodeD == 103)|(CBD_Data$DistrictCodeO == 104&CBD_Data$DistrictCodeD == 104))
DInter   = !DIntra
## Calculate the Log-likelihood function
LL = -colSums(Dtrain*log(Ptrain1) + Dbus*log(Pbus1) + Dcar*log(Pcar1) + Dbike*log(Pbike1) + Dwalk*log(Pwalk1) +
Dtrain*log(Ptrain2) + Dbus*log(Pbus2) + Dcar*log(Pcar2) + Dbike*log(Pbike2) + Dwalk*log(Pwalk2)
# DPT*log(P.PT1) + DNM*log(P.NM1) +
# DPT*log(P.PT2) + DNM*log(P.NM2) +
# DIntra*log(P.Intra) + DInter*log(P.Inter)
)
}
##### Maximize the Log-likelihood function#####
## Setting up initValaue
NL_initValue = list(Ctrain = 0,
Cbus   = 0,
Ccar   = 0,
Cbike  = 0,
TTime  = 0,
TCost  = 0,
WDayTrain = 0,
WDayBus = 0,
WDayCar = 0,
WDayBike = 0,
WDayWalk = 0,
WPurTrain = 0,
WPurBus = 0,
WPurCar = 0,
WPurBike = 0,
WPurWalk = 0,
LogsumNM = 1,
LogsumPT = 1,
LogsumInter = 1,
LogsumIntra = 1
# # Logsum1 = 1
# Logsum2 = 1
)
# Clear the Console in R
cat("\014")
## Estimate parameters
NL = mle2(NL_LL, start = NL_initValue,  control = list(trace = T, maxit = 100000000))
NL_LL0 = N*log(1/5)     ## Initial log-likelihood
NL_LL1 = logLik(NL)             ## Converged log-likelihood
## McFadden's Rho square
NL_Rho     = (1 - (NL_LL1/NL_LL0))
NL_Rho.adj = (1 - ((NL_LL1 - length(coef(NL)))/NL_LL0))
NL_GoF = cbind(NL_LL0,NL_LL1, NL_Rho, NL_Rho.adj)  ## Goodness of fit
## Preview NL model result in console
print(summary(NL))
print(NL_GoF)
library(bbmle)
source('E:/Google Drive/KVT Data/Hiroshima PhD Study/MAIN RESEARCH/2016/Summer Courses/Behavior Modeling Course 2016/BIN2016-TeamC (Hiroshima University)/BIN2016-TeamC-NL_2levels_1scaleparameter.R')
source('E:/Google Drive/KVT Data/Hiroshima PhD Study/MAIN RESEARCH/2016/Summer Courses/Behavior Modeling Course 2016/BIN2016-TeamC (Hiroshima University)/BIN2016-TeamC-NL_2levels_2scaleparameter.R')
